{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal, osr\n",
    "import os\n",
    "import io\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "from matplotlib import cm\n",
    "from matplotlib.ticker import LinearLocator, FormatStrFormatter\n",
    "import sys, os\n",
    "sys.path.append('/Volumes/EllenBackup/ANALYSIS/ANALYSIS/IS2/S2/')\n",
    "import peakdet as f1\n",
    "from scipy.stats import skew, skewtest, iqr\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from PIL import Image\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "use bash to translate TCI from jp2 to geotif (you can just run this in your terminal if you aren't using jupyter notebooks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input file size is 10980, 10980\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "!gdal_translate ~/Downloads/T10XEM_20190903T214049_TCI.jp2 ~/Downloads/T10XEM_20190903T214049_TCI.tif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "set file path (location of S2 images):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_path= '/Users/ellen.buckley/Downloads/'\n",
    "y='T10XEM_20190903T214049_TCI.tif'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "import S2 TCI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "im = gdal.Open(file_path+y)\n",
    "red=im.GetRasterBand(1).ReadAsArray()\n",
    "green=im.GetRasterBand(2).ReadAsArray() \n",
    "blue=im.GetRasterBand(3).ReadAsArray() \n",
    "\n",
    "size_m=np.shape(red)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "print info about the imported sentinel 2 image if interested:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Driver: GTiff/GeoTIFF\n",
      "Files: /Users/ellen.buckley/Downloads/T10XEM_20190903T214049_TCI.tif\n",
      "Size is 10980, 10980\n",
      "Coordinate System is:\n",
      "PROJCS[\"WGS 84 / UTM zone 10N\",\n",
      "    GEOGCS[\"WGS 84\",\n",
      "        DATUM[\"WGS_1984\",\n",
      "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n",
      "                AUTHORITY[\"EPSG\",\"7030\"]],\n",
      "            AUTHORITY[\"EPSG\",\"6326\"]],\n",
      "        PRIMEM[\"Greenwich\",0,\n",
      "            AUTHORITY[\"EPSG\",\"8901\"]],\n",
      "        UNIT[\"degree\",0.0174532925199433,\n",
      "            AUTHORITY[\"EPSG\",\"9122\"]],\n",
      "        AUTHORITY[\"EPSG\",\"4326\"]],\n",
      "    PROJECTION[\"Transverse_Mercator\"],\n",
      "    PARAMETER[\"latitude_of_origin\",0],\n",
      "    PARAMETER[\"central_meridian\",-123],\n",
      "    PARAMETER[\"scale_factor\",0.9996],\n",
      "    PARAMETER[\"false_easting\",500000],\n",
      "    PARAMETER[\"false_northing\",0],\n",
      "    UNIT[\"metre\",1,\n",
      "        AUTHORITY[\"EPSG\",\"9001\"]],\n",
      "    AXIS[\"Easting\",EAST],\n",
      "    AXIS[\"Northing\",NORTH],\n",
      "    AUTHORITY[\"EPSG\",\"32610\"]]\n",
      "Origin = (499980.000000000000000,8700000.000000000000000)\n",
      "Pixel Size = (10.000000000000000,-10.000000000000000)\n",
      "Metadata:\n",
      "  AREA_OR_POINT=Area\n",
      "Image Structure Metadata:\n",
      "  INTERLEAVE=PIXEL\n",
      "Corner Coordinates:\n",
      "Upper Left  (  499980.000, 8700000.000) (123d 0' 3.20\"W, 78d22'22.89\"N)\n",
      "Lower Left  (  499980.000, 8590200.000) (123d 0' 2.95\"W, 77d23'20.94\"N)\n",
      "Upper Right (  609780.000, 8700000.000) (118d 7'55.55\"W, 78d19'55.47\"N)\n",
      "Lower Right (  609780.000, 8590200.000) (118d30'17.01\"W, 77d21' 5.31\"N)\n",
      "Center      (  554880.000, 8645100.000) (120d39'34.59\"W, 77d52'16.59\"N)\n",
      "Band 1 Block=10980x1 Type=Byte, ColorInterp=Red\n",
      "Band 2 Block=10980x1 Type=Byte, ColorInterp=Green\n",
      "Band 3 Block=10980x1 Type=Byte, ColorInterp=Blue\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(gdal.Info(im))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "make histograms for each channel in TCI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'blue')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsYAAADSCAYAAABJsAYRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFV5JREFUeJzt3X+wXHdZx/H3h8SCQgGhUdqmNBULEjtoy6VUcRRHlLRoIw5iqoww06E6YxEGdAzClFphBnDUkaEqAUugo63FHxAhUhytg78ovZW2knaiMa0ktNgALaAIpfD4x54Lm+3ee8+9d3+cvff9mrmT3bPfu+fZk3yyzznfc3ZTVUiSJEkb3cOmXYAkSZLUBTbGkiRJEjbGkiRJEmBjLEmSJAE2xpIkSRJgYyxJkiQBNsZaQpJnJzk67TokSeqCJHclec6Q5b5frhM2xpIkSRI2xhtKks3TrkFSO+ZVkibPxnida6Z9fi3JbcD/Jnlikj9PcizJnUl+uW/sNyfZm+S+JLcDz5he5dL6lOScJB9L8oUk70nyp0levzAV2+T1U8A7m/E/nuSWJPcn+eckT+t7rlOWyPPlSa5L8u5mXQeSzE3hJUvrzTOS3N68V74zySMGBySpJN/Zd39vktf33V8015ouG+ON4SLgecDjgL8EbgVOBX4EeEWS5zbjXgc8qfl5LvDiyZcqrV9JTqCXwb308ngN8Py+IU9olp8OXJLkHOAq4BeAxwNvA/YleXiShwF/xeJ5BrgQuBZ4LLAPeOvYXpy0cfwcvffIJwFPBl67kl9eKtcjrlOrMNXGOMlVSe5N8vEWY5+Y5IbmSMttSS6YRI3rxFuq6ghwFrClqq6oqgeq6jDwdmBXM+6FwBuq6rPN+LdMqV51lJlds/OAzfQy+ZWq+gvgo32Pfw14XVV9uar+D3gp8LaqurGqvlpV7wK+3DzPM1g6zwD/WFX7q+qrwNXA94z/JaorzOvYvLWqjlTVZ4E30Dv4tBJL5VpTNu0jxnuBHS3Hvha4rqrOpvcf/++Pq6h16Ejz5+nAKc3Uzf1J7gd+Hfj25vFT+sYC/NcEa9Rs2IuZXYtTgE9WVfUt68/csar6Ut/904FXDWT2tOZ5lsszwKf6bn8ReITnLm8oezGv4zD4PnnKCn9/qVxryqbaGFfVh4HP9i9L8qQkH0xyc5J/SPJdC8OBRze3HwPcPcFSZ93Cm/AR4M6qemzfz4lVtXBk4B564VzwxIlWqc4zs2t2D3BqkvQt689cDYw/Qm8Wpz+z31JV17B8nrXBmdexGXyfHLatvgh8S9/9J/TdXirXmrJpHzEeZg/wsqp6OvArfGOv9XLgRc3nBO4HXjad8mbaR4HPNxf3fHOSTUnOSrJwkd11wKuTfGuSrbiN1Y6Zbe9fgK8ClybZnGQncO4S498O/GKSZ6bnkUmel+REls+zNIx5XbtfSrI1yePozdL86ZAxtwA/2+RyB/BDfY8tlWtNWaca4ySPAr4feE+SW+idkH5y8/BFwN6q2gpcAFzdXHyilprzDH8C+F7gTuDTwDvoHR0A+A1600J3Ah+id06itCgzuzJV9QDwU8DFwP3Ai4D30zu/cNj4eXrnI74VuA84BLykeWy5PEvHMa8j8yf03iMPNz+vHzLm5fTyeT+9i/Xeu/DAUrnW9OX4U92mUECyDXh/VZ2V5NHAwao6eci4A8CO5qIwkhwGzquqeydZr7TRmdnRSnIj8IdV9c5p16L1x7xKK9OpvcGq+jxwZ5KfBmimGBauov4EvY8jIslTgUcAx6ZSqCTAzK5Gkh9K8oTmVIoXA08DPjjturT+mVdpedP+uLZr6J1z95T0Ptj+YnpTDhcnuRU4AOxshr8KeGmz/BrgJTXtw93SBmNmR+Ip9D57+HP0ttELquqe6Zak9ci8Sis39VMpJEmSpC7o1KkUkiRJ0rTYGEuSJEn0vpp0Kk466aTatm3btFYvdc7NN9/86araMu06hjGv0vG6nFcws9KgtpmdWmO8bds25ufnp7V6qXOSdPYruM2rdLwu5xXMrDSobWY9lUKSJEnCxliSJEkCbIwlSZIkwMZYkiRJAmyMJUmSJGCKn0ohrTfbdn9g0cfueuPzJliJpDYWy6x5lbpnUu+xrY4YJ9mR5GCSQ0l2D3n8iUluSPKxJLcluWBkFUpaEfMqSdLqLNsYJ9kEXAmcD2wHLkqyfWDYa4HrqupsYBfw+6MuVNLyzKskSavX5ojxucChqjpcVQ8A1wI7B8YU8Ojm9mOAu0dXoqQVMK+SJK1Sm3OMTwWO9N0/CjxzYMzlwIeSvAx4JPCckVQnaaXMqyRJq9TmiHGGLKuB+xcBe6tqK3ABcHWShzx3kkuSzCeZP3bs2MqrlbQc8yrNGK8LkLqjTWN8FDit7/5WHjr1ejFwHUBV/QvwCOCkwSeqqj1VNVdVc1u2bFldxZKWYl6lGeJ1AVK3tGmMbwLOTHJGkhPohXLfwJhPAD8CkOSp9N5oPcQkTZ55lWaL1wVIHbJsY1xVDwKXAtcDd9Dbaz2Q5IokFzbDXgW8NMmtwDXAS6pqcPpW0piZV2nmDLsu4NSBMZcDL0pyFNgPvGzYE3n6k7R2rb7go6r20wtj/7LL+m7fDjxrtKVJWg3zKs2UlVwX8NtJvo/edQFnVdXXjvulqj3AHoC5uTl3dqVV8CuhJUmanpFdFyBp7WyMJUmaHq8LkDrExliSpCnxugCpW1qdYyxJksbD6wKk7vCIsSRJkoSNsSRJkgTYGEuSJEmAjbEkSZIE2BhLkiRJgI2xJEmSBNgYS5IkSYCNsSRJkgTYGEuSJEmAjbEkSZIE2BhLkiRJgI2xJEmSBNgYS5IkSYCNsSRJkgTYGEuSJEmAjbEkSZIE2BhLkiRJgI2xJEmSBNgYS5IkSYCNsSRJkgTYGEuSJEmAjbEkSZIE2BhLkiRJgI2xJEmSBNgYS5IkSUDLxjjJjiQHkxxKsnuRMS9McnuSA0n+ZLRlSmrLvEqStDqblxuQZBNwJfCjwFHgpiT7qur2vjFnAq8GnlVV9yX5tnEVLGlx5lWSpNVrc8T4XOBQVR2uqgeAa4GdA2NeClxZVfcBVNW9oy1TUkvmVZohzvBI3dKmMT4VONJ3/2izrN+TgScn+ackH0myY9gTJbkkyXyS+WPHjq2uYklLGVleJY1X3wzP+cB24KIk2wfG9M/wfDfwiokXKm0gbRrjDFlWA/c3A2cCzwYuAt6R5LEP+aWqPVU1V1VzW7ZsWWmtkpY3sry6IyuNnTM8Use0aYyPAqf13d8K3D1kzPuq6itVdSdwkN4br6TJGlle3ZGVxs4ZHqlj2jTGNwFnJjkjyQnALmDfwJj3Aj8MkOQkekE+PMpCJbViXqXZMbIZHnCWRxqFZRvjqnoQuBS4HrgDuK6qDiS5IsmFzbDrgc8kuR24AfjVqvrMuIqWNJx5lWbKSGdkneWR1m7Zj2sDqKr9wP6BZZf13S7glc2PpCkyr9LM+PoMD/BJejM8Pzsw5r30jhTvdYZHGj+/+U6SpClwhkfqnlZHjCVJ0ug5wyN1i0eMJUmSJGyMJUmSJMDGWJIkSQJsjCVJkiTAxliSJEkCbIwlSZIkwMZYkiRJAmyMJUmSJMDGWJIkSQJsjCVJkiTAxliSJEkCbIwlSZIkwMZYkiRJAmyMJUmSJMDGWJIkSQJsjCVJkiTAxliSJEkCbIwlSZIkwMZYkiRJAmyMJUmSJMDGWJIkSQJsjCVJkiTAxliSJEkCbIwlSZIkwMZYkiRJAmyMJUmSJMDGWJIkSQJaNsZJdiQ5mORQkt1LjHtBkkoyN7oSJa2EeZUkaXWWbYyTbAKuBM4HtgMXJdk+ZNyJwC8DN466SEntmFdp9rgzK3VHmyPG5wKHqupwVT0AXAvsHDLuN4E3A18aYX2SVsa8SjPEnVmpW9o0xqcCR/ruH22WfV2Ss4HTqur9Sz1RkkuSzCeZP3bs2IqLlbQs8yrNFndmpQ5p0xhnyLL6+oPJw4DfBV613BNV1Z6qmququS1btrSvUlJb5lWaLSPbmZW0dm0a46PAaX33twJ3990/ETgL+PskdwHnAfs8B0qaCvMqzZaR7cw6yyOtXZvG+CbgzCRnJDkB2AXsW3iwqj5XVSdV1baq2gZ8BLiwqubHUrGkpZhXabaMbGfWWR5p7ZZtjKvqQeBS4HrgDuC6qjqQ5IokF467QEntmVdp5rgzK3XI5jaDqmo/sH9g2WWLjH322suStFrmVZodVfVgkoWd2U3AVQs7s8B8Ve1b+hkkjVKrxliSJI2HO7NSd/iV0JIkSRI2xpIkSRJgYyxJkiQBNsaSJEkSYGMsSZIkATbGkiRJEmBjLEmSJAE2xpIkSRJgYyxJkiQBNsaSJEkSYGMsSZIkATbGkiRJEmBjLEmSJAE2xpIkSRJgYyxJkiQBNsaSJEkSYGMsSZIkATbGkiRJEmBjLEmSJAE2xpIkSRJgYyxJkiQBNsaSJEkSYGMsSZIkATbGkiRJEmBjLEmSJAE2xpIkSRJgYyxJkiQBLRvjJDuSHExyKMnuIY+/MsntSW5L8rdJTh99qZIkSdL4LNsYJ9kEXAmcD2wHLkqyfWDYx4C5qnoa8GfAm0ddqKR23JGVZod5lbqlzRHjc4FDVXW4qh4ArgV29g+oqhuq6ovN3Y8AW0dbpqQ23JGVZod5lbqnTWN8KnCk7/7RZtliLgb+ei1FSVo1d2Sl2WFepY5p0xhnyLIaOjB5ETAH/NYij1+SZD7J/LFjx9pXKamtke3Imldp7DzwJHVMm8b4KHBa3/2twN2Dg5I8B3gNcGFVfXnYE1XVnqqaq6q5LVu2rKZeSUsb2Y6seZXGbmR5bca4MyutUZvG+CbgzCRnJDkB2AXs6x+Q5GzgbfSa4ntHX6aklka2Iytp7EaaV3dmpbVbtjGuqgeBS4HrgTuA66rqQJIrklzYDPst4FHAe5LckmTfIk8nabzckZVmh3mVOmZzm0FVtR/YP7Dssr7bzxlxXZJWoaoeTLKwI7sJuGphRxaYr6p9HL8jC/CJqrpw0SeVNBbmVeqeVo2xpNnhjqw0O8yr1C1+JbQkSZKEjbEkSZIE2BhLkiRJgI2xJEmSBNgYS5IkSYCNsSRJkgTYGEuSJEmAjbEkSZIE2BhLkiRJgI2xJEmSBNgYS5IkSYCNsSRJkgTA5mkXMOu27f7Aoo/d9cbnTbASSeNiziVpY/CIsSRJkoRHjCVpTTyaLEnrh42xJEnSGLjjPHs8lUKSJEnCI8atLLXHJ0mStFIeTe4mG+Mx8h+9JEnS7PBUCkmSJAmPGEvS2DhrJEmzxcZYkiRplbwOaX2xMZYkSeoQZ5umx3OMJUmSJGyMJUmSJMBTKabGaRJJkqRusTGWpClw51iSusfGWJIkaUa4Uz1erRrjJDuA3wM2Ae+oqjcOPP5w4N3A04HPAD9TVXeNttSNw3/0Wgvzunp+7JKmwcxK3bFsY5xkE3Al8KPAUeCmJPuq6va+YRcD91XVdybZBbwJ+JlxFCxpceZ1fXDneOMwsxol/+9YuzZHjM8FDlXVYYAk1wI7gf7Q7gQub27/GfDWJKmqGmGtYvVHtAzEhmFe1znf+NYdM6uJ8P+Odto0xqcCR/ruHwWeudiYqnowyeeAxwOfHkWRk7Kep1HX82vTcTZMXvVQ5nwmmdkZsN6ztd5f30q0aYwzZNngXmqbMSS5BLikufs/SQ4us+6T6Ebwu1IHdKeWrtQB3all0Trypla/f/oIajCvPV2ppSt1QHdq6UodsEgtE8wrmFnoTh3QnVqs46Em8h7bpjE+CpzWd38rcPciY44m2Qw8Bvjs4BNV1R5gT5vCAJLMV9Vc2/Hj0pU6oDu1dKUO6E4tHaljw+cVulNLV+qA7tTSlTqgM7Vs+Mx2pQ7oTi3W8VCTqqXNN9/dBJyZ5IwkJwC7gH0DY/YBL25uvwD4O899kqbCvEqzxcxKHbLsEePmfKZLgevpfZTMVVV1IMkVwHxV7QP+CLg6ySF6e7G7xlm0pOHMqzRbzKzULa0+x7iq9gP7B5Zd1nf7S8BPj7Y0YAVTQmPWlTqgO7V0pQ7oTi2dqMO8At2ppSt1QHdq6Uod0JFazGxn6oDu1GIdDzWRWuJsjCRJktTuHGNJkiRp3etsY5xkR5KDSQ4l2T3hdd+V5N+S3JJkvln2uCR/k+Q/mj+/dUzrvirJvUk+3rds6LrT85ZmG92W5Jwx13F5kk822+WWJBf0Pfbqpo6DSZ47wjpOS3JDkjuSHEjy8mb5RLfJEnVMfJt0kXk1r83zdiKvy9Sy4TM7zbw2659KZruS1yVq2bCZ7VReq6pzP/QuQPhP4DuAE4Bbge0TXP9dwEkDy94M7G5u7wbeNKZ1/yBwDvDx5dYNXAD8Nb3PuDwPuHHMdVwO/MqQsdubv6OHA2c0f3ebRlTHycA5ze0TgX9v1jfRbbJEHRPfJl37Ma/mte+5O5HXZWrZ0Jmddl6bGqaS2a7kdYlaNmxmu5TXrh4x/vpXZFbVA8DCV2RO007gXc3tdwE/OY6VVNWHeejnUy627p3Au6vnI8Bjk5w8xjoWsxO4tqq+XFV3Aofo/R2Ooo57qupfm9tfAO6g9y1QE90mS9SxmLFtkw4yr+3WbV57JrVNzOxwXcwrTCCzXcnrErUsZt1ntkt57WpjPOwrMpfaQKNWwIeS3JzeNwkBfHtV3QO9v0Dg2yZYz2LrnsZ2urSZPrmqb6prInUk2QacDdzIFLfJQB0wxW3SEdN+reZ1cRs+r0NqgY2d2S68zi5ltkt5BTM79bx2tTFu9fWXY/SsqjoHOB/4pSQ/OMF1r8Skt9MfAE8Cvhe4B/jtSdWR5FHAnwOvqKrPLzV0nLUMqWNq26RDpv1azetwGz6vi9Sy0TPbhdc5C5mdxnba8JntQl672hi3+YrMsamqu5s/7wX+kt7h+f9emC5o/rx3UvUsse6Jbqeq+u+q+mpVfQ14O9+YthhrHUm+iV5Q/riq/qJZPPFtMqyOaW2TjjGvxzOvHcjrYrWY2em/zo5lthN5BTPblbx2tTFu8xWZY5HkkUlOXLgN/BjwcY7/Ss4XA++bRD2Nxda9D/j55irR84DPLUx9jMPAeUTPp7ddFurYleThSc4AzgQ+OqJ1ht63Pt1RVb/T99BEt8lidUxjm3SQeT2eeZ1yXpeqxcxOL6/Qycx2Iq+wsTPbqbzWiK/6HNUPvSsf/53elYavmeB6v4PelY63AgcW1g08Hvhb4D+aPx83pvVfQ2+64Cv09oguXmzd9KYSrmy20b8Bc2Ou4+pmPbc1/yhP7hv/mqaOg8D5I6zjB+hNj9wG3NL8XDDpbbJEHRPfJl38Ma/mtXneTuR1mVo2fGanlddm3VPLbFfyukQtGzazXcqr33wnSZIk0d1TKSRJkqSJsjGWJEmSsDGWJEmSABtjSZIkCbAxliRJkgAbY0mSJAmwMZYkSZIAG2NJkiQJgP8HUyvDfX4jZ5MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "binz=np.arange(0,270,10)\n",
    "size_m=np.shape(red)\n",
    "fig,ax=plt.subplots(1,3,figsize=(12,3))\n",
    "\n",
    "ax[0].hist(red.flatten(),bins=binz);\n",
    "ax[0].set_title('red')\n",
    "ax[1].hist(green.flatten(),bins=binz);\n",
    "ax[1].set_title('green')\n",
    "ax[2].hist(blue.flatten(),bins=binz);\n",
    "ax[2].set_title('blue')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, so clearly you can see theres a very obvious mode at bin 250-260. This is because the way the TCI is created and the bands are normalized and every value > 255 is set to 255. So these pixels are likely associated with ice (because they are so bright). So my first thought is to just set ice_pixel > 250 and see what happens:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "initiate boolean array size of S2 image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ice_mask=np.zeros(size_m).astype(bool) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "set values in the ice_mask array as True if they are greater than 250."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ice_mask[red>250]=True "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So all the remaining pixels would be open water:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "ow_mask= ~ice_mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we recreate the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_c= np.array(red,copy= True)\n",
    "green_c= np.array(green,copy= True)\n",
    "blue_c= np.array(blue,copy= True)\n",
    "\n",
    "red_c[ice_mask]= 200\n",
    "green_c[ice_mask]= 0\n",
    "blue_c[ice_mask]= 0\n",
    "\n",
    "# red_c[mp_mask]= 200\n",
    "# green_c[mp_mask]= 200\n",
    "# blue_c[mp_mask]= 0\n",
    "\n",
    "red_c[ow_mask]= 0\n",
    "green_c[ow_mask]= 0\n",
    "blue_c[ow_mask]= 200\n",
    "\n",
    "rgbArray = np.zeros((size_m[0],size_m[1],3), 'uint8')\n",
    "rgbArray[..., 0] = red_c\n",
    "rgbArray[..., 1] = green_c\n",
    "rgbArray[..., 2] = blue_c\n",
    "img = Image.fromarray(rgbArray,'RGB')\n",
    "im_save3=file_path+y[:15]+'_test_class_TCI.png'\n",
    "img.save(im_save3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That worked really well for me. the ice and open water in your image are very distinct so it was easy to separate. In sets of images with non-uniform lighting of more complex surfaces (melt ponds, ice draft, or even land...) this process gets more complicated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
